San_Juan_area
Brian
1/30/2023
Disclaimer
This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.
Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!
Methodology
SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.
Morrisey Method
The Morrisey Method is taken from Ma, Fassnacht and Kampf..
In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37
In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler
San Juan Area SNOTEL sites:
Beartown 327 Original 4/18/2005
Cascade #2 389 Morrisey 6/18/2004
Cumbres Trestle 431 Oyler -> Morrisey 6/8/2005
Idarado 538 Morrisey 7/12/2004
Lone Cone 589 Oyler -> Morrisey 6/22/2005
Middle Creek 624 Morrisey 6/26/2006
Mineral Creek 629 Oyler ?? -> Morrisey 6/22/2004
Molas Lake 632 NSCE-Morrisey, Bias-Original 10/2/2003
Red Mountain Pass 713 Morrisey 8/18/2004
Scotch Creek 739 Morrisey 6/15/2004
Slumgullion 762 Morrisey 6/26/2006
Spud Mountain 780 Morrisey 6/28/2004
Stump Lakes 797 Morrisey 7/22/2005
Upper Rio Grande 839 Oyler -> Morrisey 5/26/2004 *Why is this in red?
Upper San Juan 840 Morrisey 4/7/2004
Vallecito 843 Morrisey 7/22/2005
Wolf Creek Summit 874 Morrisey 7/12/2004
Data from thesis research:
SNOTEL_san_juan_Area <- snotel_download(site_id = c(327, 389, 431, 538, 589, 624, 629, 632, 713, 739, 762, 780, 797, 839, 840, 843, 874), path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_san_juan_Area,"C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", row.names = FALSE) #write in the raw data
Beartown 327
Original
SNOTEL_san_juan_Area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", header = TRUE)
snotel_327 <- SNOTEL_san_juan_Area %>%
filter(site_id == "327")#str(snotel_327) # check the date, usually a character.
snotel_327$Date <- as.Date(snotel_327$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_327_clean <- snotel_327 %>% # filter for the timeframe
filter(Date >= "1983-08-10" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_327_clean <- snotel_327_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_327_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_327_clean <- snotel_327_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40)ggplot(snotel_327_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_327_cull_count <- snotel_327_clean %>%
filter(temperature_min < -40) %>%
count(waterYear)
snotel_327_cull_count## # A tibble: 0 x 2
## # Groups: waterYear [0]
## # ... with 2 variables: waterYear <dbl>, n <int># filtering for too few observations in a year
snotel_327_cull_count_days <- snotel_327_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_327_cull_count_days## # A tibble: 3 x 2
## # Groups: waterYear [3]
## waterYear n
## <dbl> <int>
## 1 1983 52
## 2 1994 337
## 3 2003 3461983, 1994, 2003 need to be culled.
snotel_327_clean_culled <- snotel_327_clean %>%
filter(waterYear != "1983" & waterYear != "1994" & waterYear != "2003") #%>%
#filter(temperature_mean > -49)Beartown (327) NOT ADJUSTED.
2005-04-18 installed ext range sensor. Original
snotel_327_adjusted <- snotel_327_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-04-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
327 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_327 <- snotel_327_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_327 <- yearly_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_327 <- daily_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_327$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_327 <-daily_wy_aver_327 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_327$date_temp <- signif(daily_wy_aver2_327$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_327, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
327 SD
standard_dev_327 <- daily_wy_aver_327 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_327 <- standard_dev_327 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_327 <- standard_dev_all_327 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_327 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 8.215357 |
| 1985 | 8.122497 |
| 1986 | 7.129599 |
| 1987 | 7.619777 |
| 1988 | 8.257318 |
| 1989 | 8.044976 |
| 1990 | 7.675963 |
| 1991 | 7.684970 |
| 1992 | 7.079517 |
| 1993 | 7.642300 |
| 1995 | 7.304999 |
| 1996 | 7.681058 |
| 1997 | 7.711950 |
| 1998 | 8.105700 |
| 1999 | 6.825322 |
| 2000 | 7.526413 |
| 2001 | 7.912097 |
| 2002 | 8.099855 |
| 2004 | 7.622635 |
| 2005 | 7.225545 |
| 2006 | 7.107924 |
| 2007 | 7.426444 |
| 2008 | 7.802192 |
| 2009 | 6.893765 |
| 2010 | 7.820347 |
| 2011 | 7.632307 |
| 2012 | 7.466342 |
| 2013 | 8.000100 |
| 2014 | 7.260505 |
| 2015 | 6.671562 |
| 2016 | 7.474460 |
| 2017 | 7.072348 |
| 2018 | 6.933787 |
| 2019 | 7.658054 |
| 2020 | 7.769482 |
| 2021 | 7.841403 |
| 2022 | 7.318492 |
ggplot(standard_dev_all_327, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 327 average temperatures for water years 2005-2021
MK & SS for 327 (non-corrected)
sd_mk_327 <- mk.test(standard_dev_all_327$sd_2)
print(sd_mk_327)##
## Mann-Kendall trend test
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -140.0000000 5846.0000000 -0.2102102sd_sens_327 <- sens.slope(standard_dev_all_327$sd_2)
print(sd_sens_327)##
## Sen's slope
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.025461316 0.001275883
## sample estimates:
## Sen's slope
## -0.01223021Cascade #2 389
Morrisey 6/18/2004
snotel_389 <- SNOTEL_san_juan_Area %>%
filter(site_id == "389")#str(snotel_389) # check the date, usually a character.
snotel_389$Date <- as.Date(snotel_389$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_389_clean <- snotel_389 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_389_clean <- snotel_389_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_389_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_389_clean <- snotel_389_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_mean > -40)ggplot(snotel_389_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
temp_389_xts <- xts(snotel_389_clean$temperature_mean, order.by = snotel_389_clean$Date)
dygraph(temp_389_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_389_cull_count <- snotel_389_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_389_cull_count## # A tibble: 40 x 2
## # Groups: waterYear [40]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1984 350
## 3 1985 350
## 4 1986 365
## 5 1987 358
## 6 1988 366
## 7 1989 330
## 8 1990 348
## 9 1991 364
## 10 1992 366
## # ... with 30 more rows# filtering for too few observations in a year
snotel_389_cull_count_days <- snotel_389_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_389_cull_count_days## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1989 330
## 3 1990 348
## 4 1994 310
## 5 2006 330
## 6 2010 334
## 7 2014 314
## 8 2015 185
## 9 2016 3341983, 1989, 1990, 1994, 2006, 2010, 2014, 2015, 2016 need to be culled.
snotel_389_clean_culled <- snotel_389_clean %>%
filter(waterYear != "1983" & waterYear != "1990" & waterYear != "1994" & waterYear != "2006" & waterYear != "2010" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016") #%>%
#filter(temperature_mean > -49)Cascade #2 (389) NOT ADJUSTED.
Morrisey 6/18/2004
snotel_389_adjusted <- snotel_389_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-06-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
389 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389 <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_389 <- yearly_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_389 <- daily_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_389$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_389 <-daily_wy_aver_389 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_389$date_temp <- signif(daily_wy_aver2_389$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_389, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
389 SD
standard_dev_389 <- daily_wy_aver_389 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389 <- standard_dev_389 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389 <- standard_dev_all_389 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 10.776321 |
| 1985 | 7.879611 |
| 1986 | 10.887499 |
| 1987 | 9.408888 |
| 1988 | 10.711609 |
| 1989 | 10.418326 |
| 1991 | 10.101153 |
| 1992 | 8.637142 |
| 1993 | 8.884440 |
| 1995 | 9.371322 |
| 1996 | 9.429341 |
| 1997 | 10.027003 |
| 1998 | 9.698491 |
| 1999 | 9.029327 |
| 2000 | 9.288091 |
| 2001 | 10.630687 |
| 2002 | 10.346726 |
| 2003 | 9.727202 |
| 2004 | 9.282545 |
| 2005 | 8.472227 |
| 2007 | 9.686674 |
| 2008 | 9.541070 |
| 2009 | 8.659960 |
| 2011 | 9.334867 |
| 2012 | 9.418612 |
| 2013 | 9.751020 |
| 2017 | 9.144382 |
| 2018 | 8.911243 |
| 2019 | 9.427120 |
| 2020 | 9.880776 |
| 2021 | 9.863539 |
| 2022 | 9.441803 |
ggplot(standard_dev_all_389, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 389 average temperatures for water years 2005-2021
MK & SS for 389 (non-corrected)
sd_mk_389 <- mk.test(standard_dev_all_389$sd_2)
print(sd_mk_389)##
## Mann-Kendall trend test
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -60.0000000 3802.6666667 -0.1209677sd_sens_389 <- sens.slope(standard_dev_all_389$sd_2)
print(sd_sens_389)##
## Sen's slope
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04702986 0.01336703
## sample estimates:
## Sen's slope
## -0.01646233Corrected
389 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389_ad <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_389_ad <- yearly_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_389_ad <- daily_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_389_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_389_ad <-daily_wy_aver_389_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_389_ad$date_temp_ad <- signif(daily_wy_aver2_389_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_389_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
389 SS (corrected)
standard_dev_389_ad <- daily_wy_aver_389_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389_ad <- standard_dev_389_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389_ad <- standard_dev_all_389_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 10.061639 |
| 1985 | 7.175032 |
| 1986 | 9.930924 |
| 1987 | 8.682951 |
| 1988 | 9.868631 |
| 1989 | 9.733860 |
| 1991 | 9.472969 |
| 1992 | 8.002025 |
| 1993 | 8.283300 |
| 1995 | 8.657379 |
| 1996 | 8.697682 |
| 1997 | 9.311693 |
| 1998 | 8.929771 |
| 1999 | 8.328247 |
| 2000 | 8.510889 |
| 2001 | 9.808164 |
| 2002 | 9.563945 |
| 2003 | 8.940614 |
| 2004 | 8.587303 |
| 2005 | 8.472294 |
| 2007 | 9.686174 |
| 2008 | 9.541322 |
| 2009 | 8.659292 |
| 2011 | 9.334237 |
| 2012 | 9.418205 |
| 2013 | 9.750812 |
| 2017 | 9.144257 |
| 2018 | 8.911176 |
| 2019 | 9.426664 |
| 2020 | 9.880403 |
| 2021 | 9.863015 |
| 2022 | 9.441185 |
ggplot(standard_dev_all_389_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 389 average temperatures for water years 1986-2021
MK & SS 389 (corrected)
sd_mk_389_ad <- mk.test(standard_dev_all_389_ad$sd_2)
print(sd_mk_389_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 64.0000000 3802.6666667 0.1290323sd_sens_389_ad <- sens.slope(standard_dev_all_389_ad$sd_2)
print(sd_sens_389_ad)##
## Sen's slope
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01140359 0.04896578
## sample estimates:
## Sen's slope
## 0.01785583Cumbres Trestle 431
Oyler -> Morrisey 6/8/2005
snotel_431 <- SNOTEL_san_juan_Area %>%
filter(site_id == "431")#str(snotel_431) # check the date, usually a character.
snotel_431$Date <- as.Date(snotel_431$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_431_clean <- snotel_431 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_431_clean <- snotel_431_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_431_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_431_clean <- snotel_431_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_min < 25)ggplot(snotel_431_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_431_cull_count <- snotel_431_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_431_cull_count## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1989 364
## 4 1990 348
## 5 1991 341
## 6 1992 312
## 7 1993 356
## 8 1994 317
## 9 1995 364
## 10 1996 364
## # ... with 26 more rows# filtering for too few observations in a year
snotel_431_cull_count_days <- snotel_431_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_431_cull_count_days## # A tibble: 7 x 2
## # Groups: waterYear [7]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1990 348
## 4 1991 341
## 5 1992 312
## 6 1994 317
## 7 2014 2881987, 1988, 1990, 1991, 1992, 1994, 2014 need to be culled.
snotel_431_clean_culled <- snotel_431_clean %>%
filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_431_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_431_clean_culled <- snotel_431_clean_culled %>%
filter(temperature_mean < 20)
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Cumbres Trestle 431
Oyler -> Morrisey 6/8/2005
snotel_431_adjusted <- snotel_431_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-08", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
431 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431 <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_431 <- yearly_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_431 <- daily_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_431$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_431 <-daily_wy_aver_431 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_431$date_temp <- signif(daily_wy_aver2_431$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_431, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
431 SD
standard_dev_431 <- daily_wy_aver_431 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431 <- standard_dev_431 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431 <- standard_dev_all_431 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 8.484011 |
| 1993 | 8.101376 |
| 1995 | 7.750446 |
| 1996 | 7.952021 |
| 1997 | 8.217323 |
| 1998 | 8.600343 |
| 1999 | 7.226028 |
| 2000 | 8.164156 |
| 2001 | 8.669911 |
| 2002 | 9.005161 |
| 2003 | 8.459072 |
| 2004 | 8.306872 |
| 2005 | 8.304679 |
| 2006 | 7.417238 |
| 2007 | 7.937180 |
| 2008 | 8.259566 |
| 2009 | 7.296965 |
| 2010 | 8.418225 |
| 2011 | 8.231369 |
| 2012 | 8.081428 |
| 2013 | 8.537331 |
| 2015 | 7.219164 |
| 2016 | 7.791453 |
| 2017 | 7.322790 |
| 2018 | 7.113402 |
| 2019 | 7.880447 |
| 2020 | 8.046531 |
| 2021 | 8.095523 |
| 2022 | 7.621146 |
ggplot(standard_dev_all_431, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 431 average temperatures for water years 2005-2021
MK & SS for 431 (non-corrected)
sd_mk_431 <- mk.test(standard_dev_all_431$sd_2)
print(sd_mk_431)##
## Mann-Kendall trend test
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -104.0000000 2842.0000000 -0.2561576sd_sens_431 <- sens.slope(standard_dev_all_431$sd_2)
print(sd_sens_431)##
## Sen's slope
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.041299369 0.001003281
## sample estimates:
## Sen's slope
## -0.01890302Corrected
431 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431_ad <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_431_ad <- yearly_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_431_ad <- daily_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_431_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_431_ad <-daily_wy_aver_431_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_431_ad$date_temp_ad <- signif(daily_wy_aver2_431_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_431_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
431 SS (corrected)
standard_dev_431_ad <- daily_wy_aver_431_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431_ad <- standard_dev_431_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431_ad <- standard_dev_all_431_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 7.914872 |
| 1993 | 7.518415 |
| 1995 | 7.173944 |
| 1996 | 7.365597 |
| 1997 | 7.634730 |
| 1998 | 7.970563 |
| 1999 | 6.682639 |
| 2000 | 7.537564 |
| 2001 | 8.039524 |
| 2002 | 8.344015 |
| 2003 | 7.813684 |
| 2004 | 7.737034 |
| 2005 | 7.599212 |
| 2006 | 7.418573 |
| 2007 | 7.938448 |
| 2008 | 8.260644 |
| 2009 | 7.298209 |
| 2010 | 8.419846 |
| 2011 | 8.232506 |
| 2012 | 8.082571 |
| 2013 | 8.538590 |
| 2015 | 7.220737 |
| 2016 | 7.793079 |
| 2017 | 7.324026 |
| 2018 | 7.114853 |
| 2019 | 7.881605 |
| 2020 | 8.047598 |
| 2021 | 8.096826 |
| 2022 | 7.622159 |
ggplot(standard_dev_all_431_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 431 average temperatures for water years 1986-2021
MK & SS 431 (corrected)
sd_mk_431_ad <- mk.test(standard_dev_all_431_ad$sd_2)
print(sd_mk_431_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 54.0000000 2842.0000000 0.1330049sd_sens_431_ad <- sens.slope(standard_dev_all_431_ad$sd_2)
print(sd_sens_431_ad)##
## Sen's slope
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01334615 0.03134200
## sample estimates:
## Sen's slope
## 0.01029449Idarado 538
Morrisey 7/12/2004
snotel_538 <- SNOTEL_san_juan_Area %>%
filter(site_id == "538")#str(snotel_538) # check the date, usually a character.
snotel_538$Date <- as.Date(snotel_538$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_538_clean <- snotel_538 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_538_clean <- snotel_538_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_538_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_538_clean <- snotel_538_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)ggplot(snotel_538_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_538_cull_count <- snotel_538_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_538_cull_count## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 355
## 2 1988 366
## 3 1989 364
## 4 1990 365
## 5 1991 337
## 6 1992 358
## 7 1993 346
## 8 1994 351
## 9 1995 364
## 10 1996 365
## # ... with 26 more rows# filtering for too few observations in a year
snotel_538_cull_count_days <- snotel_538_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_538_cull_count_days## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1991 337
## 2 1993 3461991 and 1993 need to be culled.
snotel_538_clean_culled <- snotel_538_clean %>%
filter(waterYear != "1991" & waterYear != "1993") # & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_538_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_538_clean_culled <- snotel_538_clean_culled %>%
filter(temperature_mean < 20)
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Idarado 538
Morrisey 7/12/2004
snotel_538_adjusted <- snotel_538_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
538 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538 <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_538 <- yearly_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_538 <- daily_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_538$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_538 <-daily_wy_aver_538 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_538$date_temp <- signif(daily_wy_aver2_538$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_538, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
538 SD
standard_dev_538 <- daily_wy_aver_538 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538 <- standard_dev_538 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538 <- standard_dev_all_538 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.006364 |
| 1988 | 8.601864 |
| 1989 | 8.632529 |
| 1990 | 8.392986 |
| 1992 | 7.829013 |
| 1994 | 10.169815 |
| 1995 | 7.658097 |
| 1996 | 7.925011 |
| 1997 | 7.934668 |
| 1998 | 8.455869 |
| 1999 | 7.252626 |
| 2000 | 7.992580 |
| 2001 | 8.428201 |
| 2002 | 8.883209 |
| 2003 | 8.309823 |
| 2004 | 8.066471 |
| 2005 | 7.228205 |
| 2006 | 7.518896 |
| 2007 | 7.890984 |
| 2008 | 8.117130 |
| 2009 | 7.048527 |
| 2010 | 8.062733 |
| 2011 | 7.859864 |
| 2012 | 7.787906 |
| 2013 | 8.410412 |
| 2014 | 7.465408 |
| 2015 | 6.906228 |
| 2016 | 7.611126 |
| 2017 | 7.235848 |
| 2018 | 7.239780 |
| 2019 | 7.718304 |
| 2020 | 7.976082 |
| 2021 | 8.027478 |
| 2022 | 7.664387 |
ggplot(standard_dev_all_538, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 538 average temperatures for water years 2005-2021
MK & SS for 538 (non-corrected)
sd_mk_538 <- mk.test(standard_dev_all_538$sd_2)
print(sd_mk_538)##
## Mann-Kendall trend test
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -181.0000000 4550.3333333 -0.3226381sd_sens_538 <- sens.slope(standard_dev_all_538$sd_2)
print(sd_sens_538)##
## Sen's slope
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.042981556 -0.006896251
## sample estimates:
## Sen's slope
## -0.02384088Corrected
538 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538_ad <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_538_ad <- yearly_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_538_ad <- daily_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_538_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_538_ad <-daily_wy_aver_538_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_538_ad$date_temp_ad <- signif(daily_wy_aver2_538_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_538_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
538 SS (corrected)
standard_dev_538_ad <- daily_wy_aver_538_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538_ad <- standard_dev_538_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538_ad <- standard_dev_all_538_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.447584 |
| 1988 | 8.019241 |
| 1989 | 8.069873 |
| 1990 | 7.812557 |
| 1992 | 7.272544 |
| 1994 | 9.549437 |
| 1995 | 7.113641 |
| 1996 | 7.353811 |
| 1997 | 7.396458 |
| 1998 | 7.860736 |
| 1999 | 6.739299 |
| 2000 | 7.408528 |
| 2001 | 7.835100 |
| 2002 | 8.260161 |
| 2003 | 7.706375 |
| 2004 | 7.450714 |
| 2005 | 7.228187 |
| 2006 | 7.519259 |
| 2007 | 7.891254 |
| 2008 | 8.117349 |
| 2009 | 7.048722 |
| 2010 | 8.062844 |
| 2011 | 7.860215 |
| 2012 | 7.788573 |
| 2013 | 8.411054 |
| 2014 | 7.465688 |
| 2015 | 6.906947 |
| 2016 | 7.611742 |
| 2017 | 7.235817 |
| 2018 | 7.240024 |
| 2019 | 7.718445 |
| 2020 | 7.976612 |
| 2021 | 8.027693 |
| 2022 | 7.664991 |
ggplot(standard_dev_all_538_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 538 average temperatures for water years 1986-2021
MK & SS 538 (corrected)
sd_mk_538_ad <- mk.test(standard_dev_all_538_ad$sd_2)
print(sd_mk_538_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 1.100000e+01 4.550333e+03 1.960784e-02sd_sens_538_ad <- sens.slope(standard_dev_all_538_ad$sd_2)
print(sd_sens_538_ad)##
## Sen's slope
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01508416 0.01703536
## sample estimates:
## Sen's slope
## 0.0009105465Lone Cone 589
Oyler -> Morrisey 6/22/2005
snotel_589 <- SNOTEL_san_juan_Area %>%
filter(site_id == "589")#str(snotel_589) # check the date, usually a character.
snotel_589$Date <- as.Date(snotel_589$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_589_clean <- snotel_589 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_589_clean <- snotel_589_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_589_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_589_clean <- snotel_589_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)ggplot(snotel_589_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_589_cull_count <- snotel_589_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_589_cull_count## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1987 358
## 3 1988 366
## 4 1989 362
## 5 1990 365
## 6 1991 359
## 7 1992 366
## 8 1993 350
## 9 1994 337
## 10 1995 364
## # ... with 27 more rows# filtering for too few observations in a year
snotel_589_cull_count_days <- snotel_589_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_589_cull_count_days## # A tibble: 10 x 2
## # Groups: waterYear [10]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1994 337
## 3 1997 346
## 4 2006 317
## 5 2009 314
## 6 2010 332
## 7 2011 309
## 8 2012 322
## 9 2013 331
## 10 2014 3351986, 1994, 1997, 2006, 2009, 2010, 2011, 2012, 2013, 2014 need to be culled.
snotel_589_clean_culled <- snotel_589_clean %>%
filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1997" & waterYear != "2006" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_589_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_589_clean_culled <- snotel_589_clean_culled %>%
filter(temperature_mean < 20)
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Lone Cone 589
Oyler -> Morrisey 6/22/2005
snotel_589_adjusted <- snotel_589_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
589 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589 <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_589 <- yearly_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_589 <- daily_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_589$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_589 <-daily_wy_aver_589 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_589$date_temp <- signif(daily_wy_aver2_589$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_589, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
589 SD
standard_dev_589 <- daily_wy_aver_589 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589 <- standard_dev_589 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589 <- standard_dev_all_589 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.224944 |
| 1988 | 8.773281 |
| 1989 | 8.776693 |
| 1990 | 8.580289 |
| 1991 | 8.502502 |
| 1992 | 7.975611 |
| 1993 | 8.199081 |
| 1995 | 7.946318 |
| 1996 | 8.231469 |
| 1998 | 8.706109 |
| 1999 | 7.431979 |
| 2000 | 8.232988 |
| 2001 | 8.771085 |
| 2002 | 9.082003 |
| 2003 | 8.732769 |
| 2004 | 8.380419 |
| 2005 | 8.458900 |
| 2007 | 8.212611 |
| 2008 | 8.408220 |
| 2015 | 7.103666 |
| 2016 | 7.842413 |
| 2017 | 7.527960 |
| 2018 | 7.506884 |
| 2019 | 8.029199 |
| 2020 | 8.271099 |
| 2021 | 8.278625 |
| 2022 | 8.020852 |
ggplot(standard_dev_all_589, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 589 average temperatures for water years 2005-2021
MK & SS for 589 (non-corrected)
sd_mk_589 <- mk.test(standard_dev_all_589$sd_2)
print(sd_mk_589)##
## Mann-Kendall trend test
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -89.0000000 2301.0000000 -0.2535613sd_sens_589 <- sens.slope(standard_dev_all_589$sd_2)
print(sd_sens_589)##
## Sen's slope
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.043902404 0.001230005
## sample estimates:
## Sen's slope
## -0.02183399Corrected
589 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589_ad <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_589_ad <- yearly_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_589_ad <- daily_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_589_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_589_ad <-daily_wy_aver_589_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_589_ad$date_temp_ad <- signif(daily_wy_aver2_589_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_589_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
589 SS (corrected)
standard_dev_589_ad <- daily_wy_aver_589_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589_ad <- standard_dev_589_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589_ad <- standard_dev_all_589_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.671688 |
| 1988 | 8.200568 |
| 1989 | 8.222864 |
| 1990 | 7.999469 |
| 1991 | 7.962227 |
| 1992 | 7.420415 |
| 1993 | 7.639311 |
| 1995 | 7.378994 |
| 1996 | 7.641722 |
| 1998 | 8.086353 |
| 1999 | 6.914950 |
| 2000 | 7.633161 |
| 2001 | 8.156564 |
| 2002 | 8.440080 |
| 2003 | 8.103839 |
| 2004 | 7.827141 |
| 2005 | 7.781158 |
| 2007 | 8.212720 |
| 2008 | 8.408841 |
| 2015 | 7.104293 |
| 2016 | 7.842763 |
| 2017 | 7.527954 |
| 2018 | 7.507329 |
| 2019 | 8.029667 |
| 2020 | 8.271310 |
| 2021 | 8.278726 |
| 2022 | 8.021229 |
ggplot(standard_dev_all_589_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 589 average temperatures for water years 1986-2021
MK & SS 589 (corrected)
sd_mk_589_ad <- mk.test(standard_dev_all_589_ad$sd_2)
print(sd_mk_589_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.300000e+01 2.301000e+03 9.401709e-02sd_sens_589_ad <- sens.slope(standard_dev_all_589_ad$sd_2)
print(sd_sens_589_ad)##
## Sen's slope
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01153608 0.02498422
## sample estimates:
## Sen's slope
## 0.008250776Middle Creek 624
Morrisey 6/26/2006
snotel_624 <- SNOTEL_san_juan_Area %>%
filter(site_id == "624")#str(snotel_624) # check the date, usually a character.
snotel_624$Date <- as.Date(snotel_624$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_624_clean <- snotel_624 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_624_clean <- snotel_624_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_624_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_624_clean <- snotel_624_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_624_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_624_cull_count <- snotel_624_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_624_cull_count## # A tibble: 43 x 2
## # Groups: waterYear [43]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1981 354
## 3 1982 1
## 4 1983 312
## 5 1984 360
## 6 1985 362
## 7 1986 364
## 8 1987 360
## 9 1988 340
## 10 1989 363
## # ... with 33 more rows# filtering for too few observations in a year
snotel_624_cull_count_days <- snotel_624_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_624_cull_count_days## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1982 1
## 3 1983 312
## 4 1988 340
## 5 1991 338
## 6 1993 330
## 7 1994 330
## 8 2021 338
## 9 2022 3481980, 1982, 1983, 1988, 1991, 1993, 1994, 2021, 2022 need to be culled.
snotel_624_clean_culled <- snotel_624_clean %>%
filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1988" & waterYear != "1991" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022")# & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_624_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_624_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_624_clean_culled <- snotel_624_clean_culled %>%
filter(temperature_mean < 20)
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Middle Creek 624 Morrisey 6/26/2006
snotel_624_adjusted <- snotel_624_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
624 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624 <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_624 <- yearly_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_624 <- daily_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_624$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_624 <-daily_wy_aver_624 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_624$date_temp <- signif(daily_wy_aver2_624$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_624, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
624 SD
standard_dev_624 <- daily_wy_aver_624 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624 <- standard_dev_624 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624 <- standard_dev_all_624 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 8.239005 |
| 1985 | 8.207104 |
| 1986 | 7.191318 |
| 1987 | 7.871848 |
| 1989 | 8.454185 |
| 1990 | 8.200220 |
| 1992 | 6.438238 |
| 1995 | 7.565327 |
| 1996 | 7.835558 |
| 1997 | 7.996219 |
| 1998 | 8.415139 |
| 1999 | 7.174609 |
| 2000 | 7.839800 |
| 2001 | 8.355597 |
| 2002 | 8.492066 |
| 2003 | 8.150297 |
| 2004 | 7.912597 |
| 2005 | 7.915951 |
| 2006 | 7.936724 |
| 2007 | 7.640492 |
| 2008 | 8.021623 |
| 2009 | 7.171130 |
| 2010 | 8.259271 |
| 2011 | 7.932618 |
| 2012 | 7.791063 |
| 2013 | 8.259320 |
| 2014 | 7.372950 |
| 2015 | 6.806548 |
| 2016 | 7.548199 |
| 2017 | 7.108731 |
| 2018 | 7.036935 |
| 2019 | 7.675916 |
| 2020 | 7.741703 |
ggplot(standard_dev_all_624, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 624 average temperatures for water years 2005-2021
MK & SS for 624 (non-corrected)
sd_mk_624 <- mk.test(standard_dev_all_624$sd_2)
print(sd_mk_624)##
## Mann-Kendall trend test
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -128.0000000 4165.3333333 -0.2424242sd_sens_624 <- sens.slope(standard_dev_all_624$sd_2)
print(sd_sens_624)##
## Sen's slope
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0333661824 -0.0003479011
## sample estimates:
## Sen's slope
## -0.01662145Corrected
624 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624_ad <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_624_ad <- yearly_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_624_ad <- daily_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_624_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_624_ad <-daily_wy_aver_624_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_624_ad$date_temp_ad <- signif(daily_wy_aver2_624_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_624_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
624 SS (corrected)
standard_dev_624_ad <- daily_wy_aver_624_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624_ad <- standard_dev_624_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624_ad <- standard_dev_all_624_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 7.738801 |
| 1985 | 7.713059 |
| 1986 | 6.725366 |
| 1987 | 7.374899 |
| 1989 | 7.942209 |
| 1990 | 7.687289 |
| 1992 | 6.080980 |
| 1995 | 7.080976 |
| 1996 | 7.336792 |
| 1997 | 7.519955 |
| 1998 | 7.885653 |
| 1999 | 6.718142 |
| 2000 | 7.329553 |
| 2001 | 7.833709 |
| 2002 | 7.944757 |
| 2003 | 7.617514 |
| 2004 | 7.442397 |
| 2005 | 7.390605 |
| 2006 | 7.304030 |
| 2007 | 7.639661 |
| 2008 | 8.021493 |
| 2009 | 7.170919 |
| 2010 | 8.259153 |
| 2011 | 7.932368 |
| 2012 | 7.790490 |
| 2013 | 8.259010 |
| 2014 | 7.372461 |
| 2015 | 6.806595 |
| 2016 | 7.547948 |
| 2017 | 7.107845 |
| 2018 | 7.036548 |
| 2019 | 7.675791 |
| 2020 | 7.741347 |
ggplot(standard_dev_all_624_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 624 average temperatures for water years 1986-2021
MK & SS 624 (corrected)
sd_mk_624_ad <- mk.test(standard_dev_all_624_ad$sd_2)
print(sd_mk_624_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.400000e+01 4.165333e+03 4.545455e-02sd_sens_624_ad <- sens.slope(standard_dev_all_624_ad$sd_2)
print(sd_sens_624_ad)##
## Sen's slope
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01304163 0.02223675
## sample estimates:
## Sen's slope
## 0.003307866Molas Lake 632
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632 <- SNOTEL_san_juan_Area %>%
filter(site_id == "632")#str(snotel_632) # check the date, usually a character.
snotel_632$Date <- as.Date(snotel_632$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_632_clean <- snotel_632 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_632_clean <- snotel_632_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_632_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_632_clean <- snotel_632_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
#filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_632_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_632_cull_count <- snotel_632_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_632_cull_count## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1987 355
## 3 1988 362
## 4 1989 364
## 5 1990 365
## 6 1991 365
## 7 1992 366
## 8 1993 303
## 9 1994 294
## 10 1995 361
## # ... with 27 more rows# filtering for too few observations in a year
snotel_632_cull_count_days <- snotel_632_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_632_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1993 303
## 3 1994 295
## 4 2005 344
## 5 2006 304
## 6 2013 2481986, 1993, 1994, 2005, 2006, 2013 need to be culled.
snotel_632_clean_culled <- snotel_632_clean %>%
filter(waterYear != "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_632_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_632_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_632_clean_culled <- snotel_632_clean_culled %>%
filter(temperature_mean < 20)
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Molas Lake 632
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632_adjusted <- snotel_632_clean_culled %>%
mutate(temp_ad = if_else(Date < "2003-10-02", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
632 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632 <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_632 <- yearly_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_632 <- daily_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_632$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_632 <-daily_wy_aver_632 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_632$date_temp <- signif(daily_wy_aver2_632$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_632, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
632 SD
standard_dev_632 <- daily_wy_aver_632 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632 <- standard_dev_632 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632 <- standard_dev_all_632 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.613234 |
| 1988 | 8.191410 |
| 1989 | 8.166070 |
| 1990 | 7.646926 |
| 1991 | 7.969691 |
| 1992 | 7.156808 |
| 1995 | 7.489888 |
| 1996 | 7.655841 |
| 1997 | 7.698012 |
| 1998 | 8.157837 |
| 1999 | 6.885991 |
| 2000 | 7.657554 |
| 2001 | 6.788385 |
| 2002 | 8.254532 |
| 2003 | 7.907020 |
| 2004 | 7.308264 |
| 2007 | 7.607147 |
| 2008 | 7.964946 |
| 2009 | 6.957483 |
| 2010 | 7.983696 |
| 2011 | 7.701623 |
| 2012 | 7.490054 |
| 2014 | 7.236529 |
| 2015 | 6.581104 |
| 2016 | 7.409875 |
| 2017 | 7.001441 |
| 2018 | 6.814621 |
| 2019 | 7.582515 |
| 2020 | 7.766470 |
| 2021 | 7.792571 |
| 2022 | 7.357711 |
ggplot(standard_dev_all_632, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 632 average temperatures for water years 2005-2021
MK & SS for 632 (non-corrected)
sd_mk_632 <- mk.test(standard_dev_all_632$sd_2)
print(sd_mk_632)##
## Mann-Kendall trend test
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -97.0000000 3461.6666667 -0.2086022sd_sens_632 <- sens.slope(standard_dev_all_632$sd_2)
print(sd_sens_632)##
## Sen's slope
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.034372524 0.003422909
## sample estimates:
## Sen's slope
## -0.01583493Corrected
632 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632_ad <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_632_ad <- yearly_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_632_ad <- daily_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_632_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_632_ad <-daily_wy_aver_632_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_632_ad$date_temp_ad <- signif(daily_wy_aver2_632_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_632_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
632 SS (corrected)
standard_dev_632_ad <- daily_wy_aver_632_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632_ad <- standard_dev_632_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632_ad <- standard_dev_all_632_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.106416 |
| 1988 | 7.678134 |
| 1989 | 7.662133 |
| 1990 | 7.141729 |
| 1991 | 7.492618 |
| 1992 | 6.671512 |
| 1995 | 6.957618 |
| 1996 | 7.110780 |
| 1997 | 7.178235 |
| 1998 | 7.586041 |
| 1999 | 6.397699 |
| 2000 | 7.103150 |
| 2001 | 6.358660 |
| 2002 | 7.672831 |
| 2003 | 7.338249 |
| 2004 | 7.309361 |
| 2007 | 7.606845 |
| 2008 | 7.965266 |
| 2009 | 6.957636 |
| 2010 | 7.984148 |
| 2011 | 7.701741 |
| 2012 | 7.490245 |
| 2014 | 7.236308 |
| 2015 | 6.581364 |
| 2016 | 7.409996 |
| 2017 | 7.000722 |
| 2018 | 6.815246 |
| 2019 | 7.583475 |
| 2020 | 7.766635 |
| 2021 | 7.793272 |
| 2022 | 7.358050 |
ggplot(standard_dev_all_632_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 632 average temperatures for water years 1986-2021
MK & SS 632 (corrected)
sd_mk_632_ad <- mk.test(standard_dev_all_632_ad$sd_2)
print(sd_mk_632_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 57.0000000 3461.6666667 0.1225806sd_sens_632_ad <- sens.slope(standard_dev_all_632_ad$sd_2)
print(sd_sens_632_ad)##
## Sen's slope
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01086011 0.02505935
## sample estimates:
## Sen's slope
## 0.008387789